% computes Jacobian by forward differences
% Input:
% func: string, name of function
% x: point at which to take Jacobian
%    if function value at x is already known, then give x={starting point, function value}
% step: scalar, relative stepwidth;
% more arguments will be passed on to the function;
%
function jac = jacob(func,x,step,skal,varargin);
  if(iscell(x))  %function value at x is given:
    f0 = x{2};
    x = x{1};
  else
    f0 = feval(func,x,varargin{:});
  end
  n = size(x,1);
  m = size(f0,1);
  jac = zeros(m,n);
  x0 = x;
  if(isempty(skal))
    skal = ones(size(x));
  end
  for i=1:n
    step2 = step*max(1,abs(x0(i)))*skal(i);
    x = x0;
    x(i) = x0(i) + step2;
    f1 = feval(func,x,varargin{:});   
    if(any(abs(f1)>=1e100))
      error(sprintf('inadmissible value in jacob at i=%d (n=%d)',i,n));
    end
    jac(1:m,i) = (f1 - f0)/step2;
  end;
